##############################
# Analysis of survey data
##############################

rm(list=ls())

library(foreign)
library(Hmisc)
library(ri)
library(stargazer)

sink("~/Dropbox/Atacama 2015/10_replication/001_survey_analysis.txt")

##############################
# Prepare data
##############################

# load data
d = read.dta("~/Dropbox/Atacama 2015/10_replication/survey_paipote_final.dta")

# explore data
names(d)
head(d)

# gen treatment indicator
d$t_ind = NA
d$t_ind[d$a2>2]=1
d$t_ind[d$a2<3]=0
table(d$t_ind)

# gen affected area indicator 
table(d$zone)
d$zone2 = NA
d$zone2[d$zone=="A"] = 1
d$zone2[d$zone=="C"] = 1
d$zone2[d$zone=="B"] = 0
table(d$zone2)

# compare affected areas with material damage reported
table(d$t_ind,d$zone2)

# check failed conjoint experiments
table(d$failcon,d$t_ind)
reg = lm(d$failcon ~ d$t_ind)
summary(reg)

# drop fail conjoint
d = d[d$failcon==0,]

############################
# Standardized Differences
############################

# check balance for gender (placebo covariate)
round(mean(d$c1[d$t_ind==1]),2)
round(mean(d$c1[d$t_ind==0]),2)
num_c1  = abs(mean(d$c1[d$t_ind==1])- mean(d$c1[d$t_ind==0]))
den_c1 = sqrt(  ( (sd(d$c1[d$t_ind==1]))^2 + (sd(d$c1[d$t_ind==0]))^2 )/2  )
sd_c1 = abs(num_c1/den_c1)
round(sd_c1,2)

# check balance for age (placebo covariate)
round(mean(d$c2[d$t_ind==1]),2)
round(mean(d$c2[d$t_ind==0]),2)
num_c2  = abs(mean(d$c2[d$t_ind==1])- mean(d$c2[d$t_ind==0]))
den_c2 = sqrt(  ( (sd(d$c2[d$t_ind==1]))^2 + (sd(d$c2[d$t_ind==0]))^2 )/2  )
sd_c2 = abs(num_c2/den_c2)
round(sd_c2,2)

# check balance for education (placebo covariate)
round(mean(d$c3[d$t_ind==1]),2)
round(mean(d$c3[d$t_ind==0]),2)
num_c3  = abs(mean(d$c3[d$t_ind==1])- mean(d$c3[d$t_ind==0]))
den_c3 = sqrt(  ( (sd(d$c3[d$t_ind==1]))^2 + (sd(d$c3[d$t_ind==0]))^2 )/2  )
sd_c3 = abs(num_c3/den_c3)
round(sd_c3,2)

###########################
# Post treatment variables
###########################

# check ideological position
table(d$b1,exclude=NULL)
d$b1[d$b1==88]=NA
d$b1[d$b1==99]=NA
table(d$b1,exclude=NULL)
reg_ideo_con = lm(b1 ~ t_ind + c1 + c2 + c3, data = d)
summary(reg_ideo_con)

# check ideology position reported 
d$ideo_dummy = 1
d$ideo_dummy[is.na(d$b1)]=0
d$ideo_dummy[is.na(d$b1)]=0
table(d$ideo_dummy)
reg_ideo_dummy_con = lm(ideo_dummy ~ t_ind + c1 + c2 + c3, data = d)
summary(reg_ideo_dummy_con)

# Effect of floods on mayor evaluation
table(d$a3)
d$a3 = d$a3
d$a3[d$a3==88]=NA
d$a3[d$a3==99]=NA
table(d$a3)
rob5 = lm(a3 ~ t_ind + c1 + c2 + c3, data = d)
summary(rob5)

# Effect of floods on presidential evaluation
table(d$a4)
d$a4 = d$a4
d$a4[d$a4==88]=NA
d$a4[d$a4==99]=NA
table(d$a4)
rob6 = lm(a4 ~ t_ind + c1 + c2 + c3, data = d)
summary(rob6)

# check balance vote for previous president
table(d$b11)
d$bachelet = 0
d$bachelet[d$b11==1] = 1
table(d$bachelet)

round(mean(d$bachelet[d$t_ind==1]),2)
round(mean(d$bachelet[d$t_ind==0]),2)
num_bachelet  = abs(mean(d$bachelet[d$t_ind==1])- mean(d$bachelet[d$t_ind==0]))
den_bachelet = sqrt(  ( (sd(d$bachelet[d$t_ind==1]))^2 + (sd(d$bachelet[d$t_ind==0]))^2 )/2  )
sd_bachelet = abs(num_bachelet/den_bachelet)
round(sd_bachelet,2)

table(d$b11)
d$matthei = 0
d$matthei[d$b11==2] = 1
table(d$matthei)

round(mean(d$matthei[d$t_ind==1]),2)
round(mean(d$matthei[d$t_ind==0]),2)
num_matthei  = abs(mean(d$matthei[d$t_ind==1])- mean(d$matthei[d$t_ind==0]))
den_matthei = sqrt(  ( (sd(d$matthei[d$t_ind==1]))^2 + (sd(d$matthei[d$t_ind==0]))^2 )/2  )
sd_matthei = abs(num_matthei/den_matthei)
round(sd_matthei,2)

table(d$b11)
d$meo = 0
d$meo[d$b11==3] = 1
table(d$meo)

round(mean(d$meo[d$t_ind==1]),2)
round(mean(d$meo[d$t_ind==0]),2)
num_meo  = abs(mean(d$meo[d$t_ind==1])- mean(d$meo[d$t_ind==0]))
den_meo = sqrt(  ( (sd(d$meo[d$t_ind==1]))^2 + (sd(d$meo[d$t_ind==0]))^2 )/2  )
sd_meo = abs(num_meo/den_meo)
round(sd_meo,2)

table(d$b11)
d$parisi = 0
d$parisi[d$b11==4] = 1
table(d$parisi)

round(mean(d$parisi[d$t_ind==1]),2)
round(mean(d$parisi[d$t_ind==0]),2)
num_parisi  = abs(mean(d$parisi[d$t_ind==1])- mean(d$parisi[d$t_ind==0]))
den_parisi = sqrt(  ( (sd(d$parisi[d$t_ind==1]))^2 + (sd(d$parisi[d$t_ind==0]))^2 )/2  )
sd_parisi = abs(num_parisi/den_parisi)
round(sd_parisi,2)

# table ideology
stargazer(reg_ideo_con,reg_ideo_dummy_con,
          type = "text",
          title="Regression results for ideology",
          align=TRUE, 
          omit.stat=c("LL","ser","f","adj.rsq","rsq"), 
          no.space=TRUE,  
          multicolumn = TRUE,
          table.placement = "H",
          #single.row = TRUE,
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels   = c("Ideological Position", "Ideology Reported"),
          column.separate = c(1, 1),
          covariate.labels=c("Treatment"),
          dep.var.caption  = "Outcome:",
          dep.var.labels.include = FALSE,
          omit = c("c1","c2","c3","Constant"),
          add.lines = list(c("Controls","Yes","Yes")))

# table evaluations
stargazer(rob5,rob6,
          type = "text",
          title="Regression results for politicians' evaluations",
          align=TRUE, 
          omit.stat=c("LL","ser","f","adj.rsq","rsq"), 
          no.space=TRUE,  
          multicolumn = TRUE,
          table.placement = "H",
          #single.row = TRUE,
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels   = c("Mayor", "President"),
          column.separate = c(1, 1),
          covariate.labels=c("Treatment"),
          dep.var.caption  = "Performance evaluation:",
          dep.var.labels.include = FALSE,
          omit = c("c1","c2","c3","Constant"),
          add.lines = list(c("Controls","Yes","Yes")))

sink()
